↓↓ Accès au projet ↓↓
Lien GitHub du projet

Introduction

Dans le cadre du module de statistiques appliquées, nous avions abordé les bases des statistiques et pris en main des outils permettant d’effectuer des opérations statistiques avancées grâce au langage R. Ainsi, notre choix s’est porté sur des données statistiques relatives au prestige (score attribuée) de différentes professions canadiennes.

Ce projet de statistiques appliquées consiste en l’utilisation d’outils R adaptés afin d’effectuer une modélisation statistique sur notre jeu de données et produire un rapport.


Outils et environnement de travail

R et RStudio

R est un langage de programmation dont le but est de pouvoir traiter et organiser des jeux de données afin de pouvoir y appliquer des tests statistiques plus ou moins complexes et se représenter ces données graphiquement à l’aide d’une grande variété de graphiques disponibles. RStudio est une application proposant un environnement de développement et des outils adaptés au langage et à l’environnement de programmation R.

La fonction principale de RStudio consiste à faciliter le développement d’applications en langage R. Pour ce faire, le programme dispose de nombreux outils qui vous permettent notamment de créer des scripts, compiler du code, créer des graphes, ainsi que de travailler avec divers jeux de données.


R markdown

L’extension R markdown permet de générer des documents de manière dynamique en mélangeant texte mis en forme et résultats produits par du code R. Les documents générés peuvent être au format HTML, PDF, Word, et bien d’autres. C’est donc un outil très pratique pour l’exportation, la communication et la diffusion de résultats d’analyse.


Outils complémentaires

  • tidyverse : Il s’agit d’une collection d’extensions relatives à la science des données, et permettant la manipulation des tableaux de données, l’import/export de données, la manipulation de variables ou, entre autres, la visualisation de données ;
  • data.table : Extension et manipulation avancée des tableaux de données ;
  • plotly : Comme ggplot2 (compris dans le tidyverse), ce package permet la visualisation de données, à la différence près qu’il renden ces graphiques interactifs ;
  • car : Appliquer une régression linéaire et des tests statistiques ;
  • hrbrthemes : Thèmes pour ggplot2.


Données

Notre jeu de données comporte 102 individus décrits par 6 variables :


head(data) %>% paged_table()
viz_edu <- data %>% 
  ggplot(aes(x = prestige, y = education, col = type)) + 
  geom_point() + theme_bw() +
  scale_x_continuous(
    breaks = seq(25, 75, 25)
  ) +
  scale_y_continuous(
    breaks = seq(6, 16, 2)
  )

viz_edu <- ggplotly(viz_edu)

viz_inc <- data %>% 
  ggplot(aes(x = prestige, y = income, col = type)) + 
  geom_point() + theme_bw() +
  scale_x_continuous(
    breaks = seq(25, 75, 25)
  ) +
  scale_y_continuous(
    breaks = seq(0, 20000, 10000)
  )

viz_inc <- ggplotly(viz_inc)

viz_women <- data %>% 
  ggplot(aes(x = prestige, y = women, col = type)) + 
  geom_point() + theme_bw() +
  scale_x_continuous(
      breaks = seq(25, 75, 25)
    ) +
    scale_y_continuous(
      breaks = seq(0, 100, 25)
    )

viz_women <- ggplotly(viz_women)

viz_census <- data %>% 
  ggplot(aes(x = prestige, y = census, col = type)) + 
  geom_point() + theme_bw() +
  scale_x_continuous(
    breaks = seq(25, 75, 25)
  ) +
  scale_y_continuous(
    breaks = seq(2500, 7500, 2500)
  )
  
viz_census <- ggplotly(viz_census)

subplot(viz_edu, viz_inc, viz_women, viz_census, nrows = 2, margin = 0.07)

Régression linéaire

En statistiques et en apprentissage automatique, un modèle de régression linéaire est un modèle de régression qui cherche à établir une relation linéaire entre une variable, dite expliquée, et une ou plusieurs variables, dites explicatives.

Dans un premier temps, nous pouvons étudier rapidement d’éventuelles corrélations entre deux variables à l’aide d’un nuage de points. Ici, nous utilisons la fonction ggpairs afin de pouvoir en observer plusieurs en même temps.

generate_scatterM <- function(data, mapping){
  data %>% ggplot(mapping = mapping) + geom_point(size = 0.8) + 
    geom_smooth(method = loess, color = "red", size = 0.85, se = FALSE) +
    geom_smooth(method = lm, color = "blue", size = 0.85, se = FALSE)
}

scatterM <- data %>% ggpairs(columns = 1:5, lower = list(continuous = generate_scatterM))

scatterM <- ggplotly(scatterM)
scatterM

Avec les observations précédentes, nous pouvons déjà écarter certaines linéarités et choisir lesquelles sont intéressantes à étudier. Ici, nous reprenons donc la relation entre la variable prestige (il s’agit d’un score de prestige relatif à la profession) et la variable education (qui reflète le niveau d’étude). Cette fois nous utilisons ggplot pour une meilleure représentation.

education <- data %>% 
  ggplot(aes(x = education, y = prestige)) +
  geom_point() +
  geom_smooth(method = loess, se = T) +
  geom_smooth(method = lm, color = "red", se = F) +
  theme_bw() +
  scale_x_continuous(
    breaks = seq(6, 16, 2)
  ) +
  scale_y_continuous(
    breaks = seq(20, 80, 20)
  )

scatter_edu <- ggplotly(education)
scatter_edu

La droite bleue est définie par la méthode des moindres carrés (MSE), il s’agit d’une droite de régression linéaire entre les variables education et prestige.

La courbe rouge indique la tendance globale entre ces deux variables, il s’agit d’une courbe de régression de type lowess. Les deux courbes extérieures représentent l’intervalle de confiance de cette courbe de régression.

On constate que que la droite de régression est presque toujours comprise dans l’intervalle de confiance, l’hypothèse de linéarite entre les variables education et prestige est donc acceptable.


Nous testons aussi la relation entre les variables income et prestige.

income <- data %>% 
  ggplot(aes(x = income, y = prestige)) +
  geom_point() +
  geom_smooth(method = loess, se = T) +
  geom_smooth(method = lm, color = "red", se = F) +
  theme_bw() +
  scale_x_continuous(
    breaks = seq(0, 25000, 5000)
  ) +
  scale_y_continuous(
    breaks = seq(20, 100, 20)
  )

scatter_inc <- ggplotly(income)
scatter_inc

Ici, en regardant la forme du lien entre la variable entre les deux variables, on s’aperçoit que la droite de régression suis bien moins l’intervalle de confiance de la courbe lowess. l’hypothèse de linéarité est alors plus critiquable.

Pour déterminer la droite de régression, on ajuste un modèle linéaire simple aux données, à l’aide de la fonction “lm”.

lin_model_edu <- lm(prestige ~ education, data = data)
lin_model_edu
## 
## Call:
## lm(formula = prestige ~ education, data = data)
## 
## Coefficients:
## (Intercept)    education  
##     -10.732        5.361

« Intercept » correspond ici à l’ordonnée à l’origine, le « b » de notre droite, et le « x » est la pente de la droite, ce qui correspond au « b » dans notre notation.

L’équation, de notre droite est donc \(y = -10.732 + 5.361x\).


Evaluation des resultats

Le test d’évaluation de la significativité du lien linéaire entre les deux variables est valide, si les résidus :


Indépendance des résultats

On vérifie donc en premier l’hypothèse d’indépendance des résidus. On peut représenterl’autocorrélation des résidus à l’aide d’un lag_plot

bacf <- acf(residuals(lin_model_edu), plot = F)
bacfdf <- with(bacf, data.frame(lag, acf))
conf.level <- 0.95
ciline <- qnorm((1 - conf.level)/2)/sqrt(length(residuals(lin_model_edu)))

lag_plot <- bacfdf %>%
  ggplot(aes(x = lag, y = acf)) +
  geom_hline(aes(yintercept = 0)) +
  geom_segment(aes(xend = lag, yend = 0)) +
  geom_hline(aes(yintercept = ciline), linetype = 3, color = 'darkblue') +     
  geom_hline(aes(yintercept = -ciline), linetype = 3, color = 'darkblue') +
  theme_bw()

lag_plot <- ggplotly(lag_plot, tooltip = c("lag", "acf"))
lag_plot

Nous pouvons observer qu’une auto-corrélation significative est présente pour les lags 1,2,3,5,7,10 et 12. Il peut donc y avoir une auto-corrélation pour un lag de valeur 1. Nous pouvons utiliser le test de Durbin-Watson pour le confirmer.

durbinWatsonTest(lin_model_edu)
##  lag Autocorrelation D-W Statistic p-value
##    1       0.2752512       1.43917   0.004
##  Alternative hypothesis: rho != 0

Ici, le test nous indique qu’il existe une auto-corrélation significative entre les résidus d’une ligne du tableau de données et ceux de la ligne suivante. Cela signifie qu’a moins que certain residut agissent par groupe, le test rejette le lien linéaire entre les deux variables.


Distribution par loi normale


Ensuite on verifie donc que les résidut sont distribués selon une loi Normale de moyenne 0. Cette hypothèse peut s’évaluer graphiquement à l’aide d’un QQplot. Si les résidus sont bien distribués le long de la droite figurant sur le plot, alors l’hypothèse de normalité est acceptée.

qq_edu <- lin_model_edu %>%
  ggplot(aes(sample = rstandard(.))) +
  stat_qq_line(color = "red", linetype = "dashed") +
  stat_qq(size = 1.2) +
  theme_bw() +
  scale_x_continuous(
    breaks = seq(-2, 2, 1)
  ) +
  scale_y_continuous(
    breaks = seq(-3, 3, 1)
  ) +
  labs(title = "Normal Q-Q",
       x = "Theoretical Quantiles\nlm(prestige ~ education)", 
       y = "Standardized residuals")

qq_edu <- ggplotly(qq_edu)
qq_edu

Le test de Shapiro-Wilk peut également être employé pour évaluer la normalité des résidus. L’hypothèse de normalité est rejetée si la p-value est inférieure à 0.05.

shapiro.test(residuals(lin_model_edu))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(lin_model_edu)
## W = 0.98065, p-value = 0.1406


Distribution de facons homogéne


Enfin nous testons que les residus sont distribués de façon homogènes. pour cela il faut réaliser un “residuals vs fitted plot”. Les “fitted” correspondent aux réponses prédites par le modèle, pour les valeurs observées de la variable prédicitive.

sc_loc_edu <- lin_model_edu %>%
  ggplot(aes(fitted(.), sqrt(abs(rstandard(.))))) + 
  geom_point(size = 1.2) +
  geom_smooth(method = loess, se = FALSE) +
  theme_bw() +
  scale_x_continuous(
    breaks = seq(30, 70, 10)
  ) +
  scale_y_continuous(
    limits = c(0, 1.5),
    breaks = seq(0, 1.5, 0.5)
  ) +
  labs(title = "Scale-location",
       x = "Fitted values\nlm(prestige ~ education)", 
       y = "sqrt |Standardized residuals|")

sc_loc_edu <- ggplotly(sc_loc_edu)
sc_loc_edu

La courbe rouge, qui est aussi une courbe de régression locale, est globalement plate. Ceci montre que les résidus ont tendance à être répartis de façon homogène tout le long du gradient des valeurs de prestige prédites. Et donc que l’hypothèse d’homogénéité des résidus est acceptée.

Il est également possible d’évaluer cette hypothèse en employant le test de Breush-Pagan. L’hypothèse d’homogénéité est rejetée si la p-value est inférieure à 0.05.

ncvTest(lin_model_edu)
## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 0.6327545, Df = 1, p = 0.42635

Ici, le test ne rejette pas non plus l’hypothèse d’homogénéité.


Hypothése de linéarité

Il existe aussi la possibilité de directement tester la linéarité à l’aide du plot suivant.

res_v_fit_edu <- lin_model_edu %>%
  ggplot(aes(fitted(.), residuals(.))) + 
  geom_point() +
  stat_smooth(method="loess", se = FALSE) + 
  geom_hline(yintercept = 0, col = "red", linetype = "dashed") +
  theme_bw() +
  scale_x_continuous(
    breaks = seq(30, 70, 10)
  ) +
  scale_y_continuous(
    limits = c(-30, 20),
    breaks = seq(-30, 20, 10)
  ) +
  labs(title = "Residuals vs. Fitted values",
       x = "Fitted values\nlm(prestige ~ education)", 
       y = "Residuals")
    
res_v_fit_edu <- ggplotly(res_v_fit_edu)
res_v_fit_edu

Ceci nous montre que lorsque les réponses prédites par le modèle (fitted values) augmentent, les résidus restent globalement uniformément distribués de part et d’autre de 0. Cela montre, qu’en moyenne, la droite de régression, est bien adaptée aux données, et donc que l’hypothèse de linéarité est acceptable.


Interpretation des resultats

summary(lin_model_edu)
## 
## Call:
## lm(formula = prestige ~ education, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -26.0397  -6.5228   0.6611   6.7430  18.1636 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -10.732      3.677  -2.919  0.00434 ** 
## education      5.361      0.332  16.148  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.103 on 100 degrees of freedom
## Multiple R-squared:  0.7228, Adjusted R-squared:   0.72 
## F-statistic: 260.8 on 1 and 100 DF,  p-value: < 2.2e-16
confint(lin_model_edu) %>% kable()
2.5 % 97.5 %
(Intercept) -18.027220 -3.436744
education 4.702223 6.019533 
predictions <- tibble(education=c(9.21, 11.07, 14.62))
predict(lin_model_edu, newdata = predictions, interval = "confidence") %>% kable()
fit lwr upr
38.64170 36.58966 40.69374
48.61293 46.81134 50.41452
67.64405 64.52387 70.76423 
data$residuals <- residuals(lin_model_edu)
head(data) %>% kable()
education income women prestige census type residuals
13.11 12351 11.16 68.8 1113 prof 9.250875
12.26 25879 4.02 69.1 1130 prof 14.107621
12.77 9271 15.70 63.4 1171 prof 5.673573
11.42 8865 9.11 56.8 1175 prof 6.310758
14.62 8403 11.68 73.5 2111 prof 5.855950
15.64 11030 5.13 77.6 2113 prof 4.487854 
data$fitted <- fitted(lin_model_edu)
head(data) %>% kable()
education income women prestige census type residuals fitted
13.11 12351 11.16 68.8 1113 prof 9.250875 59.54913
12.26 25879 4.02 69.1 1130 prof 14.107621 54.99238
12.77 9271 15.70 63.4 1171 prof 5.673573 57.72643
11.42 8865 9.11 56.8 1175 prof 6.310758 50.48924
14.62 8403 11.68 73.5 2111 prof 5.855950 67.64405
15.64 11030 5.13 77.6 2113 prof 4.487854 73.11215 
vcov(lin_model_edu) %>% kable()
(Intercept) education
(Intercept) 13.520978 -1.1835055
education -1.183506 0.1102162 
pred_interval <- predict(lin_model_edu, interval = "prediction")
data <- cbind(data, pred_interval)
linear_regression <- data %>% 
  ggplot(aes(y = prestige, x = education)) +
  geom_point(size = 1.2) +
  geom_smooth(color = "red", method = "lm", fill = "red") +
  geom_line(aes(y = lwr), color = "red", linetype = "dashed") +
  geom_line(aes(y =  upr), color = "red", linetype = "dashed") +    
  annotate("text", x = 9, y = 80, 
           label = paste0("prestige = ", coef(lin_model_edu)["(Intercept)"] %>% round(3), 
                          " + ", coef(lin_model_edu)["education"] %>% round(3), " * education")) +
  theme_bw() +
  scale_x_continuous(
    breaks = seq(6, 16, 2)
  ) +
  scale_y_continuous(
    breaks = seq(25, 75, 25)
  ) +
  labs(title = "Linear Regression",
       x = "Education", y = "Prestige")
  

linear_regression <- ggplotly(linear_regression)
linear_regression